% =========================================================================
% Code for "Clearing Up the Fiscal Multiplier Morass"
% by Eric Leeper, Nora Traum, and Todd Walker
%
% run_prior_predictive.m
% This file lets you specificy which version of the model you would like to
% get prior-predictive government spending multipliers for, and then
% calculates the multipliers.
% =========================================================================
%
% To obtain the various scenarios displayed in Tables 3-4 and Figures 1-2
% of the main text, you must run this script for each verstion of the model
% you want to consider.

clear all; clc; close all;
rand('state',30); randn('state',30);

% -------------------------------------------------------------------------
% User Initializations 
% ------------------------------------------------------------------------- 
nprior    = 20000;              % number of draws from priors for prior predictive analysis
steps     = 1000;               % number of periods used to calculate present-value multipliers
conflevel = 0.90;               % confidence level for bands

% Select the model features you want.
% Examples for Baseline Model Scenarios:
% Model 1: All indicators = 0 
% Model 2: NK_model_I = 1; RegimeM_I = 1 for Regime M (=0 for Regime F); all other indicators zero
% Model 3: NK_model_I = 1, ROT_I = 1, RegimeM_I = 1 for Regime M (=0 for Regime F); all other indicators zero
% Model 4: NK_model_I = 1, G_in_util_I = 1 (specificy whether substitutes, complements, or unrestricted), 
%          RegimeM_I = 1 for Regime M (=0 for Regime F); all other indicators zero 
% All other variations follow from changing indicators directly
NK_model_I  = 1;               % indicator; =1 gives New Keynesian version of the model; =0 gives RBC
RegimeM_I   = 1;               % indicator; =1 gives Regime M in New Keynesian model; =0 gives Regime F
ROT_I       = 0;               % indicator; =1 gives model with nonsavers; =0 implies mu; =0 (shuts down nonsavers)
long_debt_I = 1;               % indicator; =1 gives long debt with 5-year maturity; =0 gives short debt (one period bonds)
G_in_util_I = 1;               % indicator; =1 has government spending in utility; =0 does not have G in utility function
sstax_I     = 0;               % indicator; =1 shuts down distorting tax responses; =0 implies distorting taxes may respond to debt
if G_in_util_I == 1
    Gsubstitute_I = 0;         % indicator; =1 makes public/private spending substitutes
    Gcomplement_I = 0;         % indicator; =1 makes public/private spending complements
    Gunrestricted_I = 1;       % indicator; =1 lets public/private spending be complements or substitutes
end

% -------------------------------------------------------------------------
% Run Prior Predictive
% ------------------------------------------------------------------------- 
% Define variables from DSGE model
variables;

% Initializations
PVmulti_ugc = zeros(nprior,steps,3);
IndicatorValues = [NK_model_I RegimeM_I ROT_I long_debt_I G_in_util_I sstax_I Gsubstitute_I Gcomplement_I Gunrestricted_I];
Indicators = {'NK_model_I' 'RegimeM_I' 'ROT_I' 'long_debt_I' 'G_in_util_I' 'sstax_I' 'Gsubstitute_I' 'Gcomplement_I' 'Gunrestricted_I'};

for iz = 1:nprior
    % Get prior draw of parameters
    param = get_prior_draw(IndicatorValues,Indicators);

    % Solve model for parameter draw
    [g0,g1,CC,Psi,Pie,ss,nvar,nexog] = model(param,IndicatorValues,Indicators);
    [G,C0,M,fmat,fwt,ywt,gev,eu] = gensys(g0,g1,CC,Psi,Pie);

    while eu(1,1) ~= 1 || eu(2,1) ~= 1 % Must have determinacy with a draw
        % Get prior draw of parameters
        param = get_prior_draw(IndicatorValues,Indicators);

        % Solve model for parameter draw
        [g0,g1,CC,Psi,Pie,ss,nvar,nexog] = model(param,IndicatorValues,Indicators);
        [G,C0,M,fmat,fwt,ywt,gev,eu] = gensys(g0,g1,CC,Psi,Pie);
    end

    [multi_ugc] = calculate_PVmultiplier(steps,ss,G,M,nvar,nexog);
    PVmulti_ugc(iz,:,:) = multi_ugc;
end

% -------------------------------------------------------------------------
% Calculate/Display Results 
% ------------------------------------------------------------------------- 
nperiod = 40;
X_hpd = zeros(2,nperiod,3);
Mean = zeros(nperiod,3);
for jk = 1:3 % loop through Y, C, I
    for i = 1:nperiod
        X = squeeze(PVmulti_ugc(:,i,jk));
        [X_hpd(:,i,jk)] = prob_bands(X,conflevel); 
        Mean(i,jk) = mean(squeeze(PVmulti_ugc(:,i,jk)));
    end
end
PVband = X_hpd;
PVmean = Mean;
 
Yprob = [sum(PVmulti_ugc(:,1,1) > 1)/nprior sum(PVmulti_ugc(:,4,1) > 1)/nprior sum(PVmulti_ugc(:,10,1) > 1)/nprior sum(PVmulti_ugc(:,25,1) > 1)/nprior sum(PVmulti_ugc(:,40,1) > 1)/nprior];
Cprob = [sum(PVmulti_ugc(:,1,2) > 0)/nprior sum(PVmulti_ugc(:,4,2) > 0)/nprior sum(PVmulti_ugc(:,10,2) > 0)/nprior sum(PVmulti_ugc(:,25,2) > 0)/nprior sum(PVmulti_ugc(:,40,2) > 0)/nprior];
Iprob = [sum(PVmulti_ugc(:,1,3) > 0)/nprior sum(PVmulti_ugc(:,4,3) > 0)/nprior sum(PVmulti_ugc(:,10,3) > 0)/nprior sum(PVmulti_ugc(:,25,3) > 0)/nprior sum(PVmulti_ugc(:,40,3) > 0)/nprior];

save PriorMultipliers PVmulti_ugc PVband PVmean Yprob Cprob Iprob

%% Output results in table format
FidOut = fopen('PVprobOUT.txt','w');
fprintf(FidOut,'\\begin{table}[tbp] \n');
fprintf(FidOut,'\\centering \n');
fprintf(FidOut,'\\begin{scriptsize} \n');
fprintf(FidOut,'\\begin{tabular}{llccccc} \n');
fprintf(FidOut,'\\multicolumn{7}{c}{\\textbf{Prob$\\right(PV\\frac{\\Delta Y}{\\Delta G} > 1\\left)$}}  \\\\ \n');
fprintf(FidOut,'\\noalign{\\vskip 0.1in} \n');
fprintf(FidOut,'\\hline \n');
fprintf(FidOut,'& & Impact & 4 qtrs & 10 qtrs & 25 qtrs & 10 years \\\\ \n');
fprintf(FidOut,'\\hline  \n');
fprintf(FidOut,'\\noalign{\\vskip 0.1in} \n');
fprintf(FidOut,'& & %1.2f  & %1.2f  & %1.2f  & %1.2f  & %1.2f   \\\\ \n', Yprob(1,1), Yprob(1,2), Yprob(1,3), Yprob(1,4), Yprob(1,5));
fprintf(FidOut,'& & & & & & \\\\ \n');
fprintf(FidOut,'\\multicolumn{7}{c}{\\textbf{Prob$\\right(PV\\frac{\\Delta C}{\\Delta G} > 0\\left)$}}  \\\\ \n');
fprintf(FidOut,'\\noalign{\\vskip 0.1in} \n');
fprintf(FidOut,'\\hline \n');
fprintf(FidOut,'& & Impact & 4 qtrs & 10 qtrs & 25 qtrs & 10 years \\\\ \n');
fprintf(FidOut,'\\hline  \n');
fprintf(FidOut,'\\noalign{\\vskip 0.1in} \n');
fprintf(FidOut,'& & %1.2f  & %1.2f  & %1.2f  & %1.2f  & %1.2f   \\\\ \n', Cprob(1,1), Cprob(1,2), Cprob(1,3), Cprob(1,4), Cprob(1,5));
fprintf(FidOut,'& & & & & & \\\\ \n');
fprintf(FidOut,'\\multicolumn{7}{c}{\\textbf{Prob$\\right(PV\\frac{\\Delta I}{\\Delta G} > 0\\left)$}}  \\\\ \n');
fprintf(FidOut,'\\noalign{\\vskip 0.1in} \n');
fprintf(FidOut,'\\hline \n');
fprintf(FidOut,'& & Impact & 4 qtrs & 10 qtrs & 25 qtrs & 10 years \\\\ \n');
fprintf(FidOut,'\\hline  \n');
fprintf(FidOut,'\\noalign{\\vskip 0.1in} \n');
fprintf(FidOut,'& & %1.2f  & %1.2f  & %1.2f  & %1.2f  & %1.2f   \\\\ \n', Iprob(1,1), Iprob(1,2), Iprob(1,3), Iprob(1,4), Iprob(1,5));
fprintf(FidOut,'& & & & & & \\\\ \n');
fprintf(FidOut,'\\end{tabular} \n');
fprintf(FidOut,'\\end{scriptsize} \n');
fprintf(FidOut,'\\end{table} \n');
fclose(FidOut);

%% Display multipliers
Time = 1:nperiod;
figure(1)
subplot(1,3,1)
x1 = squeeze(X_hpd(1,:,1));
x1a = squeeze(X_hpd(2,:,1));
x1b = Mean(:,1);
plot(Time,x1,'k-',Time,x1a,'k-',Time,x1b,'k-','Linewidth',1.5);
axis tight
yPos = 0.0;
hold on
plot(get(gca,'xlim'),[yPos yPos],'k:'); 
hold off
set(gca,'XLim',[0 40]);
xlabel('(a) Output Multiplier','FontSize',11)

subplot(1,3,2)
x1 = squeeze(X_hpd(1,:,2));
x1a = squeeze(X_hpd(2,:,2));
x1b = Mean(:,2);
plot(Time,x1,'k-',Time,x1a,'k-',Time,x1b,'k-','Linewidth',1.5);
axis tight
yPos = 0.0;
hold on
plot(get(gca,'xlim'),[yPos yPos],'k:'); 
hold off
set(gca,'XLim',[0 40]);
xlabel('(b) Consumption Multiplier','FontSize',11)

subplot(1,3,3)
x1 = squeeze(X_hpd(1,:,3));
x1a = squeeze(X_hpd(2,:,3));
x1b = Mean(:,3);
plot(Time,x1,'k-',Time,x1a,'k-',Time,x1b,'k-','Linewidth',1.5);
axis tight
yPos = 0.0;
hold on
plot(get(gca,'xlim'),[yPos yPos],'k:'); 
hold off
set(gca,'XLim',[0 40]);
xlabel('(c) Investment Multiplier','FontSize',11)








